An efficient scheme for numerical simulations of the spin-bath decoherence 
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We demonstrate that the Chebyshev expansion method is a very efficient numerical tool for 
studying spin-bath decoherence of quantum systems. We consider two typical problems arising 
in studying decoherence of quantum systems consisting of few coupled spins: (i) determining the 
pointer states of the system, and (ii) determining the temporal decay of quantum oscillations. As 
our results demonstrate, for determining the pointer states, the Chebyshev-based schemeis at least 
a factor of 8 faster than existing algorithms based on the Suzuki- Trotter decomposition. For the 
problems of second type, the Chebyshev-based approach has been 3-4 times faster than the Suzuki- 
Trotter-based schemes. This conclusion holds qualitatively for a wide spectrum of systems, with 
different spin baths and different Hamiltonians. 
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I. INTRODUCTION 

Recently, a great deal of attention has been devoted to 
the study of quantum computation 'l','^. For many phys- 
ical systems, basic quantum operations needed for imple- 
mentation of quantum gates have been demonstrated. To 
be practical, a quantum computer should contain a large 
number of qubits (some estimates give up to 10^ qubits 
(3), and be able to perform many hundreds of quantum 
gate operations. However, these requirements are not 
easy to satisfy in experiments. A real two-state quantum 
system is different from the ideal qubit. The system in- 
teracts with its environment, and this leads to a loss of 
phase relations between different states of the quantum 
computer (decoherence) 0, IE IE ' causing rapid accu- 
mulation of errors. Detailed theoretical understanding of 
the decoherence process is needed to prevent this. 

More generally, decoherence is an interesting many- 
body quantum phenomenon which is fundamental for 
many areas of quantum mechanics, quantum measure- 
ment theory, etc IE| - It also plays an important role in 
solid state systems, and might suppress quantum tunnel- 
ing of defects in crystals 0, spin tunneling in magnetic 
molecules and nanoparticles [E B , or destroy Kondo ef- 
fect in a dissipationless manner I.e., decoherence 
in many physical systems can have experimentally de- 
tectable (and sometimes considerable) consequences, and 
extensive theoretical studies of decoherence are needed to 
understand behavior of these systems. 

Formally speaking, decoherence is a dynamical devel- 
opment of quantum correlations (entanglement) between 
the central system and its environment. Let us assume 
that initially the central system is in the state \ipo) and 
the environment is in the state |xo)i so that the state 
of the compound system (central system plus bath) is 
\'^{t = 0)) = \ipo) (i> Ixo)- In the course of dynamical evo- 
lution, the direct product structure of the state \^{t)) 
is no longer conserved. If we need to study only the 



properties of the central system, we can consider the re- 
duced density matrix of the central system, i.e. the ma- 
trix ps{t) = Trs |*(<))(*(*)l, where Tr_B means trac- 
ing over the environmental degrees of freedom. Initially, 
Ps(0) — |V'o)('0o|; the system is in pure state, and its 
density matrix is a projector, i.e. ^1(0) = ps{0)- At 
t > 0, this property is lost, and the system appears in a 
mixed state. It has been shown that, even for relatively 
small integrable and non-integrable systems, the mixing 
is sufficient for the time-averaged, quantum dynamical 
properties of the subsystem to agree with their statistical 
mechanics values [TT|. Diagonalizing the density matrix 
ps, we can find the (instantaneous) states of the system 
\qi{t)) and (instantaneous) occupation numbers of these 
states Wi{t). It is generally assumed (and is true for all 
cases we know) that in "regular" situations, the states 
\qi{t)) quickly relax to some limiting states \pi), called 
"pointer states". This process (decoherence) is, in most 
cases, much faster than the relaxation of the occupation 
numbers (t) to their limiting values (which correspond 
to thermal equilibrium of the system with the bath). 

The theoretical description of decoherence, i.e. a de- 
scription of the evolution of the central system from its 
initial pure state ipo to the final mixed state, and finding 
the final pointer states \pi), is a very difficult problem 
of quantum many-body theory. Some simple models can 
be solved analytically, for some more complex models 
different approximations can be employed, such as the 
Markov approximation for the bath, which assumes that 
the memory effects in the bath dynamics are negligible. 
A special case of environment consisting of uncoupled 
oscillators, so-called "boson bath" , is also rather well un- 
derstood theoretically. But, although the model of boson 
bath is applicable for description of a large number of 
possible types of environments (phonons, photons, con- 
duction electrons, etc.) 0, it is not universal. 

A particularly important case where the boson bath 
description is inapplicable is the decoherence caused by 
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an enviroment made of spins, e.g. nuclear spins, or impu- 
rity spins (so called "spin bath" environment). Similarly, 
decoherence caused by some other types of quantum two- 
level systems can be described in terms of the spin bath. 
Analytical studies of the spin-bath decoherence are dif- 
ficult, and the spin-bath decoherence of many-body sys- 
tems is practically unexplored yet. In this situation, nu- 
merical modeling of spin-bath decoherence becomes an 
invaluable research tool. 

The most direct approach to study spin-bath decoher- 
ence is to compute the dynamical evolution of the whole 
compound system by directly solving the time-dependent 
Schrodinger equation of the model system. Even for a 
modest amount of spins, say 20, such calculations re- 
quire considerable computational resources, in particular 
because to study decoherence we have to follow the dy- 
namical evolution of the system over substantial periods 
of time. Therefore it is worthwhile to explore ways to 
significantly improve the efficiency of these simulations. 

In this paper, we apply the Chebyshev's expansion 
method to simulate models for the spin-bath decoher- 
ence. This method has been widely applied before 
Elllllllllllllllllll to study dynamics of large 
quantum systems, but, to our knowledge, has never been 
used for simulations of systems made of large number 
of coupled quantum spins. We show that for realistic 
problems and typical values of parameters this method 
is a very efficient tool, giving significant increase in the 
simulations speed, sometimes u p to a factor of eight, in 
comparison with the algorithms |l9ll2C| based on Suzuki- 
Trotter decompositions 21]. We illustrate this point by 
test examples that we have encountered in our previ- 
ous studies of the dynamics of the spin-bath decoher- 
ence. We also briefly discuss two other approaches, the 
short iterative Lanczos (SIL) method Clllllll and the 
multi-configurational time-dependent Hartree (MCTDH) 
[23. I25I l26i | method, which are known to demonstrate 
very good performance in many problems of quantum 
chemistry. 

The remainder of the paper is organized as follows. 
In Section we describe the model and the approaches 
used for the decoherence simulations. In Section IITTI we 
describe the specific details of application of the Cheby- 
shev's expansion method to the spin-bath decoherence 
simulations. In Section Hvl we present the results of our 
test simulations. A brief summary is given in the Section 
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II. SIMULATIONS OF THE SPIN-BATH 
DECOHERENCE: THE MODEL AND 
NUMERICAL APPROACHES 

We focus on decoherence in quantum systems of several 
coupled spins. This type of quantum systems is of partic- 
ular interest for quantum computations, since a qubit can 
be represented as a quantum spin 1/2, and qubit-based 
quantum computation is, in fact, a controlled dynamics 



of the system made of many spins 1/2. Such systems 
are also of primary interest for studying many solid state 
problems, since an electron is a particle with the spin 
1/2, and its orbital degrees of freedom are often irrele- 
vant. Thus, a system made of several coupled spins 1/2 
is a good model for investigating a large class of impor- 
tant problems both in quantum computing, and in solid 
state theory. The approach described below can be eas- 
ily extended to arbitrary spin values, but discussion of 
simulations with arbitrary spins is beyond the scope of 
this paper. 

We consider the following class of models. There is a 
central system made of AI coupled spins Sm {Sm = 1/2, 
m = 1 . . . M) . The spins Sm interact with a bath consist- 
ing of N environmental spins I„ (/„ = 1/2, n = 1 . . . N). 
The Hamiltonian governing behavior of the whole "com- 
pound" system (central spins plus the bath spins I„) 
is 

n^Ho + V = 713+^.8 + V, (1) 

where Tis and TLb are the "bare" Hamiltonians of the 
central system and the bath, correspondingly, and V is 
the system-bath interaction. Below, we present simula- 
tion results for the following general form the Hamilto- 
nians: 

■TJ \ \ ' ja QCl QCt I \ \ ^ Tja QOC 

{m,m') a=x,y,z m a=x,y,z 

{n,n') a=x,y,z n a=x,y,z 

= ^ ^ ^mn^nJn- (2) 

(m,n) a=x,y,z 

We assume that the Hamiltonian TC does not explicitly 
depend on time, i.e. all exchange interaction constants 
J, r, and A, and all external magnetic fields H are con- 
stant in time. Although this makes impossible to model 
the time-dependent quantum-gate operation, the inves- 
tigation of the fundamental properties of spin-bath de- 
coherence is not seriously affected by this requirement. 
The dynamics of the model |^ is already too complex 
to be studied analytically, and for general 7i, when no a 
priori knowledge is available, the only option is to solve 
the time-dependent Schrodinger equation of the whole 
compound system numerically. I.e., we choose some ba- 
sis states for the Hilbert space of the compound system 
(the simplest choice is the direct product of the states | t) 
and I I) for each spin S^, In)- We represent an initial 
state of the compound system 4'o as a vector in this basis 
set, and the Hamiltonian Ti. is represented as a matrix, 
so that the Schrodinger equation 

id^{t)/dt = n'^{t) (3) 

is a system of first-order ordinary differential equations 
with the initial condition '^{t = 0) = VPo- 

The length of the vector VE" is 2*^^+^; for typical values 
M = 2 and N = 20, an exact solution of about 2 • 10^ 
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differential equations becomes a serious task. Moreover, 
the interaction between the central spins is often much 
bigger than the coupling with environment or coupling 
between the bath spins, so that the system Q is often 
stiff. Simple methods, e.g. predictor-corrector schemes, 
perform rather poorly in this case, and very small inte- 
gration steps are needed to obtain a reliable solution. 

Algorithms based on the Suzuki- Trotter decomposition 
[itl i2U^ can solve (jSJ for sufhciently long times (essential 
to determine the pointer states of the central system). 
They can handle Hamiltonians with explicit dependence 
on time, are unconditionally stable, exactly preserve the 
unitarity of quantum evolution, and the time step can 
be made more than an order of magnitude bigger than 
in the typical predictor-corrector method. Moreover, as 
our experience shows, for the scheme based on Suzuki- 
Trotter decomposition, a large part of the total numerical 
error is accumulated in the total phase of the wavefunc- 
tion l^'(t)), and does not affect any measurable physical 
quantities (observables). However, for reasonably large 
systems, this scheme is still slow, and simulations of deco- 
herence lasted for up to 200 CPU hours on a SGI 3800 su- 
percomputer. The problem of long simulation times be- 
comes especially prominent if we need to find the pointer 
states, or if the dynamics of the central system is much 
faster than the decoherence rate. We found that in these 
cicumstances, the method based on Chebyshev's expan- 
sion becomes a very efficient tool to study problems of 
decoherence. 

Along with the Chebyshev's expansion method, the 
short iterative Lanczos (SIL) approach 0,|23jl23' which 
is also based on the power-series expansion of the evolu- 
tion operator, was found to be efficient for many sim- 
ilar problems of quantum chemistry. We have tested 
this method, but our results are negative. Low-order 
SIL method (with small number of Lanczos iterations 
per step, usually, less than 25) gives an unacceptable er- 
ror, even for very short time steps. On the other hand, 
high-order SIL method (with more than 25 Lanczos iter- 
ations per step) is noticeably slower the approach based 
on Chebyshev's expansion. 

We believe that low performance of SIL method orig- 
inates from the fact that for a small number of Lanc- 
zos iterations (i.e., for low-order SIL), only a very lim- 
ited part of the spectrum is described correctly. For a 
typical problem where SIL is known to be very effective 
(e.g., a wavepacket propagation), most of relevant basis 
states have energy close to the energy of a wavepacket. 
Only these relevant states should be accurately described, 
while accurate description of the whole energy spectrum 
is excessive. In contrast, in a typical spin-bath decoher- 
ence problem, a large number of bath states with very 
different energies are involved in the decoherence pro- 
cess. Correspondingly, a large part of spectrum should 
be taken into account, and the high-order SIL integrator 
should be employed, reducing the performance of the SIL 
method. 

We also note that significant speed-up can be achieved 



by using an approximate form of the wave function of 
the total system (central system plus bath). In partic- 
ular, the multi-configurational time-dependent Hartree 
(MCTDH) method Im [H |26l is known to be very effi- 
cient, e.g. for modeling of boson-bath decoherence. The 
MCTDH approach uses an approximate representation of 
the wave function, based on the assumption that the wave 
function of the total system can be written as a super- 
position of a relatively small number of " configurations" , 
i.e. products of time-varying single-spin wavefunctions. 

MCTDH is a method of choice when the dimension- 
ality of a single-particle Hilbert space is large, and the 
multi-particle quantum correlations are associated with 
a superposition of a small number of products of single- 
particle wavefunctions. The problems considered in our 
paper present an opposite situation. The bath consists of 
many spins 1/2, i.e. we have only 2 orbitals per particle 
(spin), and the single-particle evolution is very simple, 
while the complex many-particle quantum correlations 
are responsible for most of the physical effects (i.e., the 
number of important single-spin-wavefunctions products 
is very large) . It is probable that many problems of spin- 
bath decoherence can be efficiently treated by MCTDH, 
but corresponding study requires a separate extensive re- 
search effort, which is beyond the framework of our pa- 
per. 



III. CHEBYSHEV'S METHOD FOR SPIN-BATH 
DECOHERENCE 

For a time-independent Hamiltonian, the solution of 
Eq. can be formally written as 

*(t) = exp i-itn)'^Q = C/(i)*o (4) 

where U{t) = ex p j—itH ) is the evolution operator. An 
effective way [ll|ll|ll|ll|ll|ll|ll of calculation 
of the exponent of a large matrix Ti, is to expand it in 
a series of the Chebyshev polynomials of the operator 
Ti.. Below, we describe the specific details of application 
of the Chebyshev method to the spin-bath decoherence 
simulations. 

The Chebyshev's polynomials Tk{x) ~ cos (/carccosa;) 
are defined for x G [—1,1]. Thus, the Hamiltonian Ti 
first should be rescaled by the factor Eq (the range of the 
values of the system's energy) and shifted by Ec (median 
value of the systems' energy): 

Ermn = min(7^) = min ($|7i|$), 

($|$)=i 

E„,ax = niax(7^) = max ($|H|$). (5) 

($|$)=i 

In this way, the rescaled operator Q = 2{H — Ec)/Ef) 
is also bounded by —1 and 1: — 1 < (G) < 1, i.e. — 1 < 
($1^1$) < 1 for any state vector 1$) such that ($1$) = 1. 
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For spin systems, the Hamiltonian is bounded both from 
above and from below, and the operator Q can be found. 

In the specific case considered in this paper, when 
the Hamiltonian H is defined by Eq. (O, we take Eq — 
2 max {\E,nin\, \Emax\)- For this choice, ~Ea/2 < {H) < 
Eo/2. Correspondingly, we can take Ec — 0; this choice 
is legitimate, and, although might be not optimal for 
some problems, still results in very good performance of 
Chebyshev's method (see below). Since max('H) = \\H\\ 
is the norm of the Hamiltonian, the value of Eq can 
be estimated using the Cauchy's inequality: Eq/2 < 
WHsW + WUbW + M. Similarly, 

W-HsW < E E l-^w|-|l^™M|5™'ll (6) 

{m,m') a=x,y,z 

+E E 

m a—x.y,z 

""EE 4i^wi+E E 2''^™'' 

{711, m') Oi—x.y,z m a—x.y,z 

and llTisll and || V|| can be estimated in the same manner. 
As a result, we have an estimate Eq < Ei, where 

■^1 = E E 2'^^™™''^ 51 2'^""'' 

{m.m') a=x,y,z {n,7i') 

+ E ^K»J+Ei^™i+Ei^"i' (7) 

and the operator Q can be defined a,s G — 27Y / Ei , which 
satisfies the inequality — 1 < (G) < 1- 

The Chebyshev's expansion of the evolution operator 
U{t) (see Eq. 0J now looks like 

CxD 

U{t) - exp (~^rg) = E ""^Tk {G) (8) 

fc=0 

where r = Eit/2. The expansion coefficients Ck can be 
calculated using the orthogonal property of the polyno- 
mials Tk{x): 

ak f'^ Tk{x)exp{-ixT) •\fe 7 / \ fn\ 

Ck^ — / , dx = ak[-i) Jk{T), (9) 

TT J_i VI - x^ 

where Jfc(T) is the Bessel function of /c-th order, and ak — 
2 for fc = and Ofc = 1 for k >1. The successive terms 
in the Chebyshev's series can be efficiently determined 
using the recursion 

Tk+i{G) = 2GTk{G) + Tk-i{G) (10) 

with the conditions To{G) = 1, Ti{G) ^ G- Thus, to 
find the vector ^'(i), we just need to sum successively 
the terms of the series (O, using Eq. (|10() for calculation 
of the subsequent terms, until we reach some pre-defined 
value K oi k, which is determined by the required preci- 
sion. 




FIG. 1: Dependence of the order of the Chebyshev's expan- 
sion K on the value of r = Eit/2. The solid circles cor- 
responds to the minimum value e = IG^"' of the expansion 
coefficient Ck\ the open circles corresponds to e = 10~^. 

The high precision of this scheme originates from the 
fact that, for k ^ t, the value of a Bessel function de- 
creases super-exponentially Jfc(T) ~ {r/k)^ , so that ter- 
mination of the series ^ sX k = K leads to an error 
which decreases super-exponentially with K. In prac- 
tice, K = 1.5t already gives precision of 10~^ or bet- 
ter in most cases. Due to the same reason, this scheme 
is asymptotically more efficient than any time-marching 
scheme. For given sufficiently small error e, the num- 
ber of operations Nop needed for finding the wavefunc- 
tion at time T, i.e. '^{T), grows linearly with T for 
the Chebyshev-based scheme. For a marching scheme 
of order r with the time step At, the numerical error 
is e ~ {AtyX, so that for given e and T, the number 
of operations needed is Nop = T/At ^ j^i+i/r^ grow- 
ing super-linear ly with increasing T. For very long-time 
simulations, and when very high precision is necessary, 
the Chebyshev method is more efficient than any time- 
marching scheme known to us. However, in practice, 
a precision better than 0.5%~1% is very rarely needed. 
Similarly, very long-time simulations are rarely of inter- 
est; in most cases, the simulations are interesting only 
until the dynamics of the system exhibits some non- 
trivial behavior. Therefore, in spite of its asymptotic 
efficiency, the Chebyshev method is not always the best 
choice for real research, and its efficiency should be stud- 
ied in every separate case. 



IV. SIMULATION RESULTS 

We assess the usefulness of the Chebyshev's method for 
a wide spectrum of decoherence problems, by consideri- 
ong two central problems of decoherence, description of 
damping of quantum oscillations in a system, and deter- 
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mination of the pointer states. In fact, there is no strict 
boundary: studying both problems, we track evolution of 
the system checking its state at regular intervals of length 
T, but in studying the oscillations decay the interval T 
is much smaller than the characteristic decoherence time 
Tdec, while in studying the pointer states, T is larger than 

In spite of the asymptotic advantages of the 
Chebyshev-based scheme, it is not a priori clear if it is 
efficient for realistic problems, when the required numer- 
ical error 5 is modest (say, 5 = 10~^-10~'^). Also, if we 
track the dynamics of the decoherence process, we make 
many steps of modest length T, and the overhead asso- 
ciated with the use of the Chebyshev's expansion might 
be significant, see Fig.Q] 

To study this issue, we have performed several types 
of numerical tests. The timing information reported in 
this paper has been obtained from calculations on a SGI 
Origin 3800 (500 MHz) system, using sequential, single 
processor code. The order of Chebyshev's expansion K 
have been defined by the pre-specified precision e. We 
determined the minimum value of K such that \ck\ < e 
for k > K, starting from the value Ko ~ [I-It] {[x] is 
the integer part of x), and adjusting it as needed. Each 
simulation has been performed three times: (i) using the 
Chebyshev's method with e = 10^^^, the reference run, 

(ii) using Chebyshev's method with e ~ 10~^-10~^, and 

(iii) usin g t he scheme based on Suzuki- Trotter decom- 
position [igl . Previously we have used the latter to 
study spin-bath decoherence ^10, 27] . In this paper, we 
have chosen to consider the same problems as in our pre- 
vious works on this subject, in order to avoid the impres- 
sion that the tests have been constructed to favor one 
particular method. 

First, we consider the problem of oscillations damping 
in the central system of two spins coupled by Heisenberg 
exchange, interacting with the bath. We studied this 
problem using the Suzuki- Trotter scheme in Ref.l2^ The 
Hamiltonians describing the bath and the system are: 

Hs = JSiS2, Hb^O, V = ^A„(Si + S2)I„, 

(11) 

with = 16 bath spins. The exchange parameter 
J — 16.0 (antiferromagnetic coupling between the cen- 
tral spins), while An are uniformly distributed between 
and —0.5. The initial state of the compound system 
1^*0) = IV'o) ^ Ixo) is the product of the initial state \ipo) 
of the central system, and |xo) of the bath. In this case, 
l''/'o) = I Ti)? i-e. the first central spin is in the state 
S^{t = 0) = +1/2, and the second spin is in the state 
Sf{t = 0) = -1/2. The initial state of the bath |xo) is 
the linear superposition of all basis states with random 
coefficients. Physically, this situation corresponds to the 
case of the temperature 6 which is high in comparison 
with the bath energies An, but is much lower than the 
system's energy J (note that J ^ An in this case). 

The initial state of the central system is a superpo- 
sition of two eigenstates of H: the state with the total 
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FIG. 2: Time dependence of the oscillations of the expectation 
value of Sl{t) in the two-spin system decohered by a spin 
bath. 



spin S = 1 and Sz = 0, and the state with the total 
spin S = 0. These states have different energies, and, for 
example, the dynamics of (t) is represented by oscilla- 
tions with the frequency J. Due to interaction with the 
spin bath, these oscillations are damped, see Fig. [3 To 
study this damping in detail, we take the Suzuki- Trotter 
time step At = 0.035, T — 2At, and watch the system 
since t = till t^ax = SOOT. If we do not need such a 
high resolution, we increase T. In Table we present 
the CPU time needed to perform the simulations using 
the Suzuki- Trotter and Chebyshev's methods, along with 
the resulting error S (which should not be confused with 
the "nominal" precision of the Chebyshev's scheme e). 
The error S has been obtained from comparison with the 
"reference" Chebyshev's run (e = 10""'^^), and is equal to 
the maximum of absolute errors of the quantities (all nor- 
malized to unity) 2S'f, 2S^, 45"? S'fja, /3 = x,y,z), and 
the so-called "quadratic entropy" [23 5*'^' = 1 — Trp|. 
These quantities have been calculated and compared at 
regular intervals of length T. Their calculation increases 
the number of computations, so that the tests 1, 2, and 
3, which are otherwise equivalent for the Suzuki- Trotter 
method, require more and more CPU time. 



TABLE I: Comparison of the Suzuki- Trotter scheme (abbre- 
viated as ST) with the Chebyshev's scheme (abbreviated as 
Ch) for the problem of oscillations decay. 



Test 


At 


T 


tmax 




S 




e 


CPU Time 


1, Ch 




200Af 


ST 


1 • 


IQ-^ 




10-*^ 


22 min 


1, ST 


0.035 


200Af 


8T 


0.44 


■ 10" 


-2 




SO min 


2, Ch 




SAt 


200T 


0.3 


10" 


4 


lO"*^ 


59 min 


2, ST 


0.035 


8At 


20or 


0.48 


• 10" 


-2 




89 min 


3, Ch 




2At 


SOOT 


0.55 


■ 10" 


-3 


10"'^ 


226 min 


3, ST 


0.035 


2At 


SOOT 


0.48 


■ 10" 


-2 




156 min 
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TABLE II: Comparison of the Suzuki- Trotter scheme (abbre- 
viated as ST) with the Chebyshev's scheme (abbreviated as 
Ch) for the problem of oscillations decay, employing the two- 
leap approach with different Ti and T2. 



Test 


At 


Ti 


T-z 




S 




e 


CPU Time 


4, Ch 




150Af 


At 


0.4 


■ IQ- 


'4 


10"^ 


61 min 


4, ST 


0.02 


150Af 


At 


0.2 


■ 10" 


-2 




144 min 


5, Ch 




300Af 


At 


0.4 


• 10" 


-4 


lO"*^ 


75 min 


5, ST 


0.02 


SOOAf 


At 


0.3 


■ 10" 


-2 




221 min 



As one can see from Table 13 for realistic values of 
maximum error S ~ 0.5 • 10~^, and even for not very 
long runs, the Chebyshev's scheme can be faster than 
the Suzuki- Trotter method by a factor of up to four, and 
the efficiency of the Chebyshev's scheme grows fast with 
increasing T. However, this straightforward comparison 
is too crude, and Table is only an illustration of basic 
features of the Chebyshev's method. To model fast oscil- 
lations which decay slowly (often, with the decay time of 
order of decoherence time Tdec), we should make T signif- 
icantly smaller than the oscillation period tosc = 27r/ J, in 
order to correctly determine the amplitude of oscillations 
at given time. 

Therefore, to track the damping of oscillations, we use 
the two- leap approach: first, we make a large time leap 
of length Ti (Ti > tosc, but Ti < Tdec), and then we 
make a number n2 (usually, 15-20) of smaller steps T2 
such that T2 <ti tosc but n2T2 > tosc, resolving in detail 
one period of oscillations and extracting the amplitude. 
By repeating this two-stage sequence ntot times, we can 
reliably track the change of the oscillations amplitude 
with time. The test example of this approach have been 
taken from our recent work [2^ . We have performed the 
same kind of simulations as described above, with N — 
16 bath spins, repeating the two-leap sequence ntot — 8 
times, each time making one long leap Ti followed by 
7T.2 = 21 short leaps T2. The results of these tests are 
presented in Table IIT| Again, Chebyshev-based method 
can be up to three times faster than the Suzuki- Trotter 
algorithm [lll20l|. 

Finally, we have tested the Chebyshev scheme in the 
problem of determining the pointer states, employing an 
example from our work IQJ. This example is interesting 
also because it deals with a physically important case of 
a spin bath possessing chaotic internal dynamics, which 
is relevant for majority of realistic spin baths (such as 
nuclear spins or impurity spins baths). The Hamiltonian 
describing the system is 



(12) 



i.e., the bath spins are coupled only with the first central 
spin, and the bath Hamiltonian is now 
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500 1000 1500 2000 

Time 

FIG. 3: Temporal evolution of different elements of the den- 
sity matrix p: diagonal elements corresponding to the states 
TT), I Ti), I it), and I li) (the four upper curves), and the 
non-diagonal element p]^ (the lowest curve). Very slow re- 
laxation is better seen for the uppermost curve (the diagonal 
element corresponding to the state | ||)) which has a small 
negative slope. Note that the two lines in the middle (the 
second and the third lines from above, the diagonal elements 
corresponding to the states | jj,) and | It), correspondingly) 
are very close to each other at t > 200, as expected for a near- 
equilibrium (although not completely relaxed) situation. 



In our simulations we used h = 0.1 and £/„„/ randomly 
distributed in the interval [-0.013,0.013]. This Hamil- 
tonian is known to result in stochastic behavior |30j : we 
have checked the level statistics independently, and found 
that it closely follows the Wigner-Dyson distribution. 

To determine the pointer states, we need to find the 
elements of the reduced density matrix ps (t) in the long- 
time limit t ^ 00. We start at t = from the state 
of the compound system which is the product of the 
states of the bath and the central system (as above), 
but the initial state of the central spins now is the sin- 
glet li^o) - (l/%/2)[| Ti) - I it)]- Because of decoher- 
ence, the final state of the central system is mixed, and 
Ps = wi\pi){pi\ + 'W2\p2) {p2\, where \pi) and \p2) are 
the pointer states, which are superpositions of the states 
I TT)jJii) I Ti)) and I iT). As we have found in our 
work |lOj| . the form of this superposition is determined 
by the ratio J/b, where b = '^n^n- For J/b ^ 1, the 
pointer states are very close to the singlet 5 = and 
triplet S — I, Sz — states, and for J <C &, the pointer 
states are close to | Ti) and | iT). Thus, the quantities 
characterizing the type of the pointer state are the values 
of the non-diagonal elements of the density matrix ps in 
the basis | TT), I ii) I Ti), and | iT). In particular, the 
element = (Ti \ps\ iT) is a very suitable quantity 
to characterize the pointer state. This non-diagonal ele- 
ment is close to zero for J <C 6, and gradually increases 
in absolute value with increasing J. 

Typical results for temporal evolution of the elements 
of the density matrix ps are shown in Fig. |3| One can 
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TABLE III: Comparison of the Suzuki- Trotter scheme (ab- 
breviated as ST) with the Chebyshev's scheme (abbreviated 
as Ch) for the problem of determining the pointer states. 
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see that in this situation, we do not need to use the two- 
leap approach with different Ti and T2. The relaxation 
(after some initial period) is slow, and no fast oscilla- 
tions of considerable amplitude exist at long times, so 
that the one-leap approach is sufficient. Thus, the ef- 
ficiency of the Chebyshev-based scheme is expected to 
be very good. This is indeed the Table UTTI 

demonstrates. The results presented there correspond 
to J = 0.1. The Chebyshev-based scheme is faster than 
the Suzuki- Trotter method up to a factor of 8. 

We have checked our conclusions on many other cases, 
with the central systems made of up to M = 4 spins, and 
with the baths made of up to = 22 spins, with different 
Hamiltonians and different values of the Hamiltonian pa- 
rameters. We found that Chebyshev-based method gives 
a significant increase in the simulations speed for all prob- 
lems where the value of T can be made sufficiently large. 

V. SUMMARY 

Theoretical studies of the spin-bath decoherence are 
important for many areas of physics, including quantum 



mechanics and quantum measurement theory, quantum 
computing, sohd state physics etc. Decoherence is a com- 
plex many-body phenomenon, and numerical simulation 
is an important tool for its investigation. In this paper, 
we have studied efficiency of the numerical scheme based 
on the Chebyshev expansion. We have presented specific 
details of the application of this method to the spin-bath 
decoherence modeling. To assess the efficiency of the 
simulation method, we have used model problems which 
we have encountered in our previous studies of the spin- 
bath decoherence. We compared the Chebyshev-based 
scheme with a fast method based on the Suzuki- Trotter 
decomposition. We have found that in many cases, the 
former gives a considerable increase in the speed of simu- 
lations, sometimes up to a factor of eight (for the problem 
of finding the system's pointer states), while in studying 
the decoherence dynamics, the increase in speed is less 
drastic (a factor of 2-3), but still considerable. This con- 
clusion holds for many types of central systems and spin 
baths, with different Hamiltonians. 
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